##########################################################################################################################################
###### This file estimates and saves all the betas used
###############################

# set to your working directory:
setwd("C:\\Users\\Admin\\OneDrive\\Desktop\\Market Power and Systematic Risk - Data and Code")
getwd()

#To build the VAR, we need the factors described in Campbell and Vuolteenaho (2004):

#four state variables: excess log market return, term yield spread, Shiller's Price Earnings (PE) and VS

library(zoo)
library(stringi)
library(data.table)
## excess log market return:##

######
# these files are downloaded from Kenneth French's data library
VWret = read.table("VWMktRet.txt", sep="\t",dec=".",header=T)
#Combine with FF, to get the risk-free rates:
FFRF = read.table("C:F-F_Research_Data_Factors.txt",sep="",dec=".",header=T)

VWret$DATE = as.yearmon(as.character(VWret$DATE),format="%Y%m%d")
FFRF$Date = as.yearmon(as.character(FFRF$Date),format="%Y%m")

VWret = merge(VWret,FFRF,by.x="DATE",by.y="Date")

#######
# This file is downloaded from Amit Goyal's webpage
TY = read.table("PredictorData2019.csv",dec=",",header=T,sep=";")

TY = TY[,c("yyyymm","lty","tbl","Rfree")]
#In annualized percentage points
TY$TY = TY$lty - TY$tbl
TY = na.omit(TY)
TY$TY = TY$TY*100
TY$yyyymm = as.yearmon(as.character(TY$yyyymm),format="%Y%m")

#######
# This file is downloaded from Robert Shiller's webpage
#PE - ratio
PE = read.table("ie_data.csv",sep=";",dec=",",header=T)
PE = PE[,c(1,13)]
PE$Date = as.numeric(paste0(format(seq(from=as.Date("1871-01-01",format="%Y-%d-%m"),to=as.Date("2020-01-06",format="%Y-%d-%m"),by="month"),"%Y"),
                            format(seq(from=as.Date("1871-01-01",format="%Y-%d-%m"),to=as.Date("2020-01-06",format="%Y-%d-%m"),by="month"),"%m")))

PE$Date = as.yearmon(as.character(PE$Date),format="%Y%m")

PE$CAPE = log(PE$CAPE)

#####################################
# this is a file containing data from the CRSP/Compustat merged database. It includes
# the following data columns:
# Date, PERMNO, permco, cfacpr, cfacshr, shrout, altprc, prc, ret, retx, vol,
# shrcd, exchcd, siccd, dlstcd, dlret, retadj.1mn, ME, first, last, port.weight,
# datadate, comp.count, at, revt, ib, sich, dvc, BE, OpProf, GrProf, Cflow, Inv,
# AstChg, Davis.bkeq
# Note that some of these elements are modified raw data, calculating the variables
# of Fama and French (see the description in the French Data Library)
load("200426 data.both.m.RData")
#########################################

#Historical BE's in Davis.bkeq

#Use the Breakpoints:
######
# these files are from Kenneth French's data library 
#First column tells us how many are below zero, then the number of observations left
Break_BEME = read.table("BE-ME_Breakpoints.txt",sep="",header=F)
#First column tells us how many are below zero, then the number of observations left
Break_ME = read.table("ME_Breakpoints.txt",sep="",header=F)
#One Breakpoint for Break_ME (the median -- > which is column 12)
#Two Breakpoints for BEME, 30% and 70%, which is column 8 and 16
Break_BEME = Break_BEME[,c(1,8,16)]

#Pick out all small stocks with high BE/ME annually (end of June)
#picks all size related content at June
Break_ME = Break_ME[,c(1,12)][stri_paste_list(lapply(strsplit(as.character(Break_ME[,c(1)]),split=""),function(x) tail(x,n=2, collapse = ''))) == "06",]
#Delete the 2020 observation, because it is not available in BEME
Break_ME = Break_ME[-nrow(Break_ME),]

#Convert to Day-Count Convention:
Break_ME$V1 = seq(1926.417,2019.417,by=1)
Break_BEME$V1 = seq(1926.417,2019.417,by=1)

#now these are the correct breakpoints, now select BE and ME from a DEC in t-1

VS = data.frame(Date = unique(data.both.m$Date), VS = NA)
#Calculate the small-stock value spread each June:
for(i in 1:nrow(Break_BEME)){
  #Select the matching dates:
  Size_Low = data.both.m[data.both.m$Date == as.yearmon(Break_ME$V1[i]),]$ME<Break_ME$V12[i]
  BE_Interm = c()
  #Select which companies are in which portfolio:
  BE_Interm = data.both.m[data.both.m$Date == as.yearmon(Break_BEME$V1[i]),]$BE  
  BE_Interm[is.na(BE_Interm)] = data.both.m[data.both.m$Date == as.yearmon(Break_BEME$V1[i]),]$Davis.bkeq[is.na(BE_Interm)]
  BEME_Low = (BE_Interm/data.both.m[data.both.m$Date == as.yearmon(Break_BEME$V1[i]),]$ME)<Break_BEME$V8[i]
  
  BEME_High = (BE_Interm/data.both.m[data.both.m$Date == as.yearmon(Break_BEME$V1[i]),]$ME)>Break_BEME$V16[i]
  
  #Difference between the Intersection
  VS$VS[VS$Date == as.yearmon(Break_BEME$V1[i])] = mean(log(BE_Interm[BEME_High & Size_Low& BE_Interm > 0]/data.both.m[data.both.m$Date == as.yearmon(Break_BEME$V1[i]),]$ME[BEME_High & Size_Low& BE_Interm > 0]),na.rm=T) -
    mean(log(BE_Interm[BEME_Low & Size_Low & BE_Interm > 0]/data.both.m[data.both.m$Date == as.yearmon(Break_BEME$V1[i]),]$ME[BEME_Low & Size_Low& BE_Interm > 0]),na.rm=T)
}


######
# these files are from Kenneth French's data library 
#Load the FF 2x3 Portfolio:
Six_Portfolios = read.table("6_Portfolios_2x3.txt",header=T,sep="")
#Transform the equally weighted portfolio into log returns:
Transform_Simple = function(x){
  Log_Ret = c()
  Log_Ret = log(1+x)
  return(Log_Ret)
}

Six_Portfolios[,2:length(Six_Portfolios)] = apply(Six_Portfolios[,2:length(Six_Portfolios)]/100,2,Transform_Simple)

Six_Portfolios$Date = as.yearmon(as.character(Six_Portfolios$Date), format="%Y%m")
Filled_Months = which(!is.na(VS$VS))
for(i in 1:(length(Filled_Months)-1)){
  #From the first one to the last one
  Temp = VS$VS[Filled_Months[i]:Filled_Months[i+1]]
  Temp_Dates = VS$Date[Filled_Months[i]:Filled_Months[i+1]]
  Temp_Dates = Temp_Dates[-c(1,13)]
  SLoBM = Six_Portfolios$SmallLoBM[Six_Portfolios$Date%in%Temp_Dates]
  SHiBM = Six_Portfolios$SmallHiBM[Six_Portfolios$Date%in%Temp_Dates]
  
  Temp[2:(length(Temp)-1)] = Temp[1] + cumsum(SLoBM - SHiBM)
  VS$VS[Filled_Months[i]:Filled_Months[i+1]] = Temp
}

#Merge: VS VWret, TY, PE:
#in order: R_Me, TY, PE, VS
CheckCheck = merge(VWret,TY,by.x="DATE",by.y="yyyymm")
CheckCheck = merge(CheckCheck,PE,by.x="DATE",by.y="Date")
CheckCheck = merge(CheckCheck,VS,by.x="DATE",by.y="Date")
colnames(CheckCheck)[c(2,11,12,13)] = c("R_Me","TY","PE","VS")

RF = CheckCheck$Rfree

CheckCheck$R_Me = CheckCheck$R_Me - CheckCheck$Rfree
CheckCheck = CheckCheck[,c("DATE","R_Me","TY","PE","VS")]

#Gamma: A Matrix with all coefficient estimates, the top row are the estimates from the return regressions, in Campbell rho:

Gamma_Check = matrix(NA,nrow=4,ncol=4)
#Recreate the VAR(1) estimates --> works:
Gamma_Check[1,] = coef(lm(CheckCheck[,2]~shift(CheckCheck[,2],n=1)+shift(CheckCheck[,3],n=1)+shift(CheckCheck[,4],n=1)+shift(CheckCheck[,5],n=1)))[2:5]

Gamma_Check[2,] = coef(lm(CheckCheck[,3]~shift(CheckCheck[,2],n=1)+shift(CheckCheck[,3],n=1)+shift(CheckCheck[,4],n=1)+shift(CheckCheck[,5],n=1)))[2:5]

Gamma_Check[3,] = coef(lm(CheckCheck[,4]~shift(CheckCheck[,2],n=1)+shift(CheckCheck[,3],n=1)+shift(CheckCheck[,4],n=1)+shift(CheckCheck[,5],n=1)))[2:5]

Gamma_Check[4,] = coef(lm(CheckCheck[,5]~shift(CheckCheck[,2],n=1)+shift(CheckCheck[,3],n=1)+shift(CheckCheck[,4],n=1)+shift(CheckCheck[,5],n=1)))[2:5]


#Residual Matrix:
Residual_Check = matrix(NA,ncol=(nrow(CheckCheck)-1),nrow=4)

Residual_Check[1,] = c(residuals(lm(CheckCheck[,2]~shift(CheckCheck[,2],n=1)+shift(CheckCheck[,3],n=1)+shift(CheckCheck[,4],n=1)+shift(CheckCheck[,5],n=1))))

Residual_Check[2,] = c(residuals(lm(CheckCheck[,3]~shift(CheckCheck[,2],n=1)+shift(CheckCheck[,3],n=1)+shift(CheckCheck[,4],n=1)+shift(CheckCheck[,5],n=1))))

Residual_Check[3,] = c(residuals(lm(CheckCheck[,4]~shift(CheckCheck[,2],n=1)+shift(CheckCheck[,3],n=1)+shift(CheckCheck[,4],n=1)+shift(CheckCheck[,5],n=1))))

Residual_Check[4,] = c(residuals(lm(CheckCheck[,5]~shift(CheckCheck[,2],n=1)+shift(CheckCheck[,3],n=1)+shift(CheckCheck[,4],n=1)+shift(CheckCheck[,5],n=1))))

rho = 0.95^(1/12) #rho is set to 0.95!! annually (...)

#e1 is a vector with a 1 in the first element and zero anywhere else:

Market.Forecast = c(NA,(lm(CheckCheck[,2]~shift(CheckCheck[,2],n=1)+shift(CheckCheck[,3],n=1)+shift(CheckCheck[,4],n=1)+shift(CheckCheck[,5],n=1)))$fitted)

Identity_Matrix = matrix(0, nrow=4, ncol=4)
diag(Identity_Matrix) = rep(1,4)

Lambda_Check = rho * Gamma_Check %*% solve(Identity_Matrix - rho*Gamma_Check)
#e1
e1 = c(1, 0, 0, 0)

N_CF = (e1 + e1%*%Lambda_Check)%*%(Residual_Check)
N_DR = e1%*%Lambda_Check%*%Residual_Check

#Add an NA to N_CF and N_DR

N_CF = c(NA,N_CF)
N_DR = c(NA,-N_DR)


library(rmatio)
###############
# This is a file containing the monthly returns from CRSP, adjusted for delistings
# following Warther (1995) and Shumway and Warther (1999)
# the file format is the following:
# first column: date identifiers (yyyymmdd)
# first row: PERMNO identifiers
# elements 2:end,2:end: the returns in decimals
# First part of the list: monthly CRSP value weighted return
# Second part of the list: monthly S&P 500 return
# Third part of the list: Returns of individual stocks
CrossSection=read.mat("MR_ab1926(dlr).mat")
# this file has a similar setup. It contains the stock price and market capitalization 
# information from CRSP
# First part of the list: stock prices
# Second part of the list: market capitalization
Value = read.mat("Pr_MK.mat")
###############

#Remove all PERMNO that are in CrossSection, but not in Value!
CrossSection[[3]] = CrossSection[[3]][,CrossSection[[3]][1,]%in%Value[[2]][1,]]
Sec_ID = CrossSection[[3]][1,]
Value[[2]] = Value[[2]][,Value[[2]][1,]%in%CrossSection[[3]][1,]]

CrossSection[[3]] = CrossSection[[3]][CrossSection[[3]][,1]<20190101,]
CrossSection[[3]] = CrossSection[[3]][CrossSection[[3]][,1]>19260701,]

CrossSection[[2]] = CrossSection[[2]][CrossSection[[2]][,1]<20190101,]
CrossSection[[2]] = CrossSection[[2]][CrossSection[[2]][,1]>19260701,]

CrossSection[[1]] = CrossSection[[1]][CrossSection[[1]][,1]<2019101,]
CrossSection[[1]] = CrossSection[[1]][CrossSection[[1]][,1]>19260701,]

CrossSection[[3]] = CrossSection[[3]][-1,]


#Idiosyncratic Volatility

###
# this file is from Kenneth French's data library
FF = read.table("F-F_Research_Data_Factors.txt",sep="",header=T)
#Reduce the Matrix (...)
FF = FF[FF$Date<201901,]
FF = FF[FF$Date>=192607,]
FF$Mkt.RF = FF$Mkt.RF/100
FF$SMB = FF$SMB/100
FF$HML = FF$HML/100
FF$RF = FF$RF/100


##################################### Efficient Solution ###################################################

library(data.table)
library(roll)
#They all have the same length already (have to do nothing!)


#Calculate excess returns
Monthly_Data = apply(CrossSection[[3]],2,function(x) x-RF)
#convert this to a dataframe
Monthly_Data = data.frame(Monthly_Data)
#merge all variables on to it (Mkt, N_DR, N_CF)
Mkt_Monthly = CheckCheck$R_Me #and N_DR, N_CF
Monthly_Data$R_Me = Mkt_Monthly
Monthly_Data$N_DR = N_DR
Monthly_Data$N_CF = N_CF
#convert it to a data table
setDT(Monthly_Data)
#Traditional Beta Estimation (use a 5-year, 60 months window:)
#Market
cols <- colnames(Monthly_Data)[2:(ncol(Monthly_Data)-3)]
Beta_M = matrix(NA,nrow=dim(Monthly_Data)[1],ncol=(dim(Monthly_Data)[2]-3))
Beta_M[,2:(dim(Monthly_Data)[2]-3)] = as.matrix(Monthly_Data[ , lapply(.SD, function(X) roll_cov(cbind(X,R_Me),width = 60,min_obs = 24)[1,2,]/roll_cov(cbind(R_Me,R_Me),width = 60,min_obs = 24)[1,1,]),.SDcols=cols])
#Up
Monthly_Data_Up = as.data.frame(Monthly_Data)
setDT(Monthly_Data_Up)
v1 <- colnames(Monthly_Data_Up)[2:(ncol(Monthly_Data_Up)-3)]
for(j in v1){
  set(Monthly_Data_Up, i = which(Monthly_Data_Up$R_Me < 0), j = j, value = NA_integer_)
}
Beta_Up = matrix(NA,nrow=dim(Monthly_Data)[1],ncol=(dim(Monthly_Data)[2]-3))
Beta_Up[,2:(dim(Monthly_Data_Up)[2]-3)] = as.matrix(Monthly_Data_Up[ , lapply(.SD, function(X) roll_cov(cbind(X,R_Me),width = 60,min_obs = 12)[1,2,]/roll_cov(cbind(R_Me,R_Me),width = 60,min_obs = 12)[1,1,]),.SDcols=cols])

#Down
Monthly_Data_Down = as.data.frame(Monthly_Data)
setDT(Monthly_Data_Down)
v1 <- colnames(Monthly_Data_Down)[2:(ncol(Monthly_Data_Down)-3)]
for(j in v1){
  set(Monthly_Data_Down, i = which(Monthly_Data_Down$R_Me > 0), j = j, value = NA_integer_)
}
Beta_Down = matrix(NA,nrow=dim(Monthly_Data)[1],ncol=(dim(Monthly_Data)[2]-3))
Beta_Down[,2:(dim(Monthly_Data_Down)[2]-3)] = as.matrix(Monthly_Data_Down[ , lapply(.SD, function(X) roll_cov(cbind(X,R_Me),width = 60,min_obs = 12)[1,2,]/roll_cov(cbind(R_Me,R_Me),width = 60,min_obs = 12)[1,1,]),.SDcols=cols])

#CF
cols <- colnames(Monthly_Data)[2:(ncol(Monthly_Data)-3)]
Beta_CF = matrix(NA,nrow=dim(Monthly_Data)[1],ncol=(dim(Monthly_Data)[2]-3))
Beta_CF[,2:(dim(Monthly_Data)[2]-3)] = as.matrix(Monthly_Data[ , lapply(.SD, function(X) roll_cov(cbind(X,N_CF),width = 60,min_obs = 24)[1,2,]/roll_cov(cbind(N_CF-N_DR,N_CF-N_DR),width = 60,min_obs = 24)[1,1,]),.SDcols=cols])
#DR
cols <- colnames(Monthly_Data)[2:(ncol(Monthly_Data)-3)]
Beta_DR = matrix(NA,nrow=dim(Monthly_Data)[1],ncol=(dim(Monthly_Data)[2]-3))
Beta_DR[,2:(dim(Monthly_Data)[2]-3)] = as.matrix(Monthly_Data[ , lapply(.SD, function(X) roll_cov(cbind(X,N_DR),width = 60,min_obs = 24)[1,2,]/roll_cov(cbind(N_CF-N_DR,N_CF-N_DR),width = 60,min_obs = 24)[1,1,]),.SDcols=cols])

#Other Beta Estimation Methods:
#1. Weighting:
h = log(2)/(2/3*60)
Weights = exp(-abs(60:1)*h)/(sum(exp(-abs(60:1)*h)))

#Market
cols <- colnames(Monthly_Data)[2:(ncol(Monthly_Data)-3)]
Weighting_Beta_M = matrix(NA,nrow=dim(Monthly_Data)[1],ncol=(dim(Monthly_Data)[2]-3))
Weighting_Beta_M[,2:(dim(Monthly_Data)[2]-3)] = as.matrix(Monthly_Data[ , lapply(.SD, function(X) roll_cov(cbind(X,R_Me),width = 60,min_obs = 24,weights=Weights)[1,2,]/roll_cov(cbind(R_Me,R_Me),width = 60,min_obs = 24,weights=Weights)[1,1,]),.SDcols=cols])
#Up
Monthly_Data_Up = as.data.frame(Monthly_Data)
setDT(Monthly_Data_Up)
v1 <- colnames(Monthly_Data_Up)[2:(ncol(Monthly_Data_Up)-3)]
for(j in v1){
  set(Monthly_Data_Up, i = which(Monthly_Data_Up$R_Me < 0), j = j, value = NA_integer_)
}
Weighting_Beta_Up = matrix(NA,nrow=dim(Monthly_Data)[1],ncol=(dim(Monthly_Data)[2]-3))
Weighting_Beta_Up[,2:(dim(Monthly_Data_Up)[2]-3)] = as.matrix(Monthly_Data_Up[ , lapply(.SD, function(X) roll_cov(cbind(X,R_Me),width = 60,min_obs = 12,weights=Weights)[1,2,]/roll_cov(cbind(R_Me,R_Me),width = 60,min_obs = 12,weights=Weights)[1,1,]),.SDcols=cols])

#Down
Monthly_Data_Down = as.data.frame(Monthly_Data)
setDT(Monthly_Data_Down)
v1 <- colnames(Monthly_Data_Down)[2:(ncol(Monthly_Data_Down)-3)]
for(j in v1){
  set(Monthly_Data_Down, i = which(Monthly_Data_Down$R_Me > 0), j = j, value = NA_integer_)
}
Weighting_Beta_Down = matrix(NA,nrow=dim(Monthly_Data)[1],ncol=(dim(Monthly_Data)[2]-3))
Weighting_Beta_Down[,2:(dim(Monthly_Data_Down)[2]-3)] = as.matrix(Monthly_Data_Down[ , lapply(.SD, function(X) roll_cov(cbind(X,R_Me),width = 60,min_obs = 12,weights=Weights)[1,2,]/roll_cov(cbind(R_Me,R_Me),width = 60,min_obs = 24,weights=Weights)[1,1,]),.SDcols=cols])

#CF
cols <- colnames(Monthly_Data)[2:(ncol(Monthly_Data)-3)]
Weighting_Beta_CF = matrix(NA,nrow=dim(Monthly_Data)[1],ncol=(dim(Monthly_Data)[2]-3))
Weighting_Beta_CF[,2:(dim(Monthly_Data)[2]-3)] = as.matrix(Monthly_Data[ , lapply(.SD, function(X) roll_cov(cbind(X,N_CF),width = 60,min_obs = 24,weights=Weights)[1,2,]/roll_cov(cbind(N_CF-N_DR,N_CF-N_DR),width = 60,min_obs = 24,weights=Weights)[1,1,]),.SDcols=cols])
#DR
cols <- colnames(Monthly_Data)[2:(ncol(Monthly_Data)-3)]
Weighting_Beta_DR = matrix(NA,nrow=dim(Monthly_Data)[1],ncol=(dim(Monthly_Data)[2]-3))
Weighting_Beta_DR[,2:(dim(Monthly_Data)[2]-3)] = as.matrix(Monthly_Data[ , lapply(.SD, function(X) roll_cov(cbind(X,N_DR),width = 60,min_obs = 24,weights=Weights)[1,2,]/roll_cov(cbind(N_CF-N_DR,N_CF-N_DR),width = 60,min_obs = 24,weights=Weights)[1,1,]),.SDcols=cols])

#save the final dataframes already, so we can use them for the Shrinkage:

library(data.table)
library(lubridate)


Beta_M[,1] = Monthly_Data[[1]]
#Remove the last three columns, where we have saved R_ME, N_DR and N_CF before:
Beta_M = data.frame(Beta_M)
library(stringr)
colnames(Beta_M)[1:ncol(Beta_M)] = c("Date",str_remove_all(cols,pattern="X"))
Beta_M$Date = as.Date(as.character(Beta_M$Date),format="%Y%m%d")


Beta_Up[,1] = Monthly_Data[[1]]
#Remove the last column, where we have saved vwret before:
Beta_Up = data.frame(Beta_Up)
colnames(Beta_Up)[1:ncol(Beta_Up)] = c("Date",str_remove_all(cols,pattern="X"))
Beta_Up$Date = as.Date(as.character(Beta_Up$Date),format="%Y%m%d")


Beta_Down[,1] = Monthly_Data[[1]]
#Remove the last column, where we have saved vwret before:
Beta_Down = data.frame(Beta_Down)
colnames(Beta_Down)[1:ncol(Beta_Down)] = c("Date",str_remove_all(cols,pattern="X"))
Beta_Down$Date = as.Date(as.character(Beta_Down$Date),format="%Y%m%d")


Beta_CF[,1] = Monthly_Data[[1]]
#Remove the last column, where we have saved vwret before:
Beta_CF = data.frame(Beta_CF)
colnames(Beta_CF)[1:ncol(Beta_CF)] = c("Date",str_remove_all(cols,pattern="X"))
Beta_CF$Date = as.Date(as.character(Beta_CF$Date),format="%Y%m%d")


Beta_DR[,1] = Monthly_Data[[1]]
#Remove the last column, where we have saved vwret before:
Beta_DR = data.frame(Beta_DR)
colnames(Beta_DR)[1:ncol(Beta_DR)] = c("Date",str_remove_all(cols,pattern="X"))
Beta_DR$Date = as.Date(as.character(Beta_DR$Date),format="%Y%m%d")


Weighting_Beta_M[,1] = Monthly_Data[[1]]
#Remove the last column, where we have saved vwret before:
Weighting_Beta_M = data.frame(Weighting_Beta_M)
colnames(Weighting_Beta_M)[1:ncol(Weighting_Beta_M)] = c("Date",str_remove_all(cols,pattern="X"))
Weighting_Beta_M$Date = as.Date(as.character(Weighting_Beta_M$Date),format="%Y%m%d")


Weighting_Beta_Up[,1] = Monthly_Data[[1]]
#Remove the last column, where we have saved vwret before:
Weighting_Beta_Up = data.frame(Weighting_Beta_Up)
colnames(Weighting_Beta_Up)[1:ncol(Weighting_Beta_Up)] = c("Date",str_remove_all(cols,pattern="X"))
Weighting_Beta_Up$Date = as.Date(as.character(Weighting_Beta_Up$Date),format="%Y%m%d")


Weighting_Beta_Down[,1] = Monthly_Data[[1]]
#Remove the last column, where we have saved vwret before:
Weighting_Beta_Down = data.frame(Weighting_Beta_Down)
colnames(Weighting_Beta_Down)[1:ncol(Weighting_Beta_Down)] = c("Date",str_remove_all(cols,pattern="X"))
Weighting_Beta_Down$Date = as.Date(as.character(Weighting_Beta_Down$Date),format="%Y%m%d")


Weighting_Beta_CF[,1] = Monthly_Data[[1]]
#Remove the last column, where we have saved vwret before:
Weighting_Beta_CF = data.frame(Weighting_Beta_CF)
colnames(Weighting_Beta_CF)[1:ncol(Weighting_Beta_CF)] = c("Date",str_remove_all(cols,pattern="X"))
Weighting_Beta_CF$Date = as.Date(as.character(Weighting_Beta_CF$Date),format="%Y%m%d")


Weighting_Beta_DR[,1] = Monthly_Data[[1]]
#Remove the last column, where we have saved vwret before:
Weighting_Beta_DR = data.frame(Weighting_Beta_DR)
colnames(Weighting_Beta_DR)[1:ncol(Weighting_Beta_DR)] = c("Date",str_remove_all(cols,pattern="X"))
Weighting_Beta_DR$Date = as.Date(as.character(Weighting_Beta_DR$Date),format="%Y%m%d")

#save:
save(Beta_M,file="Beta_M.RData")
save(Beta_CF,file="Beta_CF.RData")
save(Beta_DR,file="Beta_DR.RData")
save(Beta_Up,file="Beta_Up.RData")
save(Beta_Down,file="Beta_Down.RData")
save(Weighting_Beta_M,file="Weight_Beta_M.RData")
save(Weighting_Beta_CF,file="Weight_Beta_CF.RData")
save(Weighting_Beta_DR,file="Weight_Beta_DR.RData")
save(Weighting_Beta_Up,file="Weight_Beta_Up.RData")
save(Weighting_Beta_Down,file="Weight_Beta_Down.RData")

#2. Shrinkage:

load("Beta_M.RData")
load("Beta_Up.RData")
load("Beta_Down.RData")
load("Beta_CF.RData")
load("Beta_DR.RData")
#This is the Cross-section
Shrinking_Beta_M = matrix(NA,nrow=dim(Beta_M)[1],ncol=dim(Beta_M)[2])
Shrinking_Beta_Up = matrix(NA,nrow=dim(Beta_M)[1],ncol=dim(Beta_M)[2])
Shrinking_Beta_Down = matrix(NA,nrow=dim(Beta_M)[1],ncol=dim(Beta_M)[2])
Shrinking_Beta_CF = matrix(NA,nrow=dim(Beta_M)[1],ncol=dim(Beta_M)[2])
Shrinking_Beta_DR = matrix(NA,nrow=dim(Beta_M)[1],ncol=dim(Beta_M)[2])
library(data.table)
setDF(Beta_M)
setDF(Beta_Up)
setDF(Beta_Down)
setDF(Beta_CF)
setDF(Beta_DR)
#For Beta_M
cross_Beta_M = rowMeans(Beta_M[,2:ncol(Beta_M)], na.rm=T)
sd_cross_Beta_M = (sd(rowMeans(Beta_M[,2:ncol(Beta_M)], na.rm=T),na.rm=T))^2
prior_Beta_M = (apply(Beta_M[,2:ncol(Beta_M)],1,sd,na.rm=T))^2
#For Beta_Up
cross_Beta_Up = rowMeans(Beta_Up[,2:ncol(Beta_Up)], na.rm=T)
sd_cross_Beta_Up = (sd(rowMeans(Beta_Up[,2:ncol(Beta_Up)], na.rm=T),na.rm=T))^2
prior_Beta_Up = (apply(Beta_Up[,2:ncol(Beta_Up)],1,sd,na.rm=T))^2
#For Beta_Down
cross_Beta_Down = rowMeans(Beta_Down[,2:ncol(Beta_Down)], na.rm=T)
sd_cross_Beta_Down = (sd(rowMeans(Beta_Down[,2:ncol(Beta_Down)], na.rm=T),na.rm=T))^2
prior_Beta_Down = (apply(Beta_Down[,2:ncol(Beta_Down)],1,sd,na.rm=T))^2
#For Beta CF
cross_Beta_CF = rowMeans(Beta_CF[,2:ncol(Beta_CF)], na.rm=T)
sd_cross_Beta_CF = (sd(rowMeans(Beta_CF[,2:ncol(Beta_CF)], na.rm=T),na.rm=T))^2
prior_Beta_CF = (apply(Beta_CF[,2:ncol(Beta_CF)],1,sd,na.rm=T))^2
#For Beta DR
cross_Beta_DR = rowMeans(Beta_DR[,2:ncol(Beta_DR)], na.rm=T)
sd_cross_Beta_DR = (sd(rowMeans(Beta_DR[,2:ncol(Beta_DR)], na.rm=T),na.rm=T))^2
prior_Beta_DR = (apply(Beta_DR[,2:ncol(Beta_DR)],1,sd,na.rm=T))^2

for(i in 2:ncol(Beta_M)){
  Shrinking_Beta_M[,i] = (Beta_M[,i] + sd_cross_Beta_M/prior_Beta_M * cross_Beta_M)/(1 + sd_cross_Beta_M/prior_Beta_M)
  
  Shrinking_Beta_Up[,i] = (Beta_Up[,i] + sd_cross_Beta_Up/prior_Beta_Up * cross_Beta_Up)/(1 + sd_cross_Beta_Up/prior_Beta_Up)
  
  Shrinking_Beta_Down[,i] = (Beta_Down[,i] + sd_cross_Beta_Down/prior_Beta_Down * cross_Beta_Down)/(1 + sd_cross_Beta_Down/prior_Beta_Down)
  
  Shrinking_Beta_CF[,i] = (Beta_CF[,i] + sd_cross_Beta_CF/prior_Beta_CF * cross_Beta_CF)/(1 + sd_cross_Beta_CF/prior_Beta_CF)
  
  Shrinking_Beta_DR[,i] = (Beta_DR[,i] + sd_cross_Beta_DR/prior_Beta_DR * cross_Beta_DR)/(1 + sd_cross_Beta_DR/prior_Beta_DR)
  
  print(paste0(i," of iteration ",ncol(Beta_M)," which is ",i/ncol(Beta_M)*100," Percent"))
}
#Add the Dates:
Shrinking_Beta_M = data.frame(Shrinking_Beta_M)
Shrinking_Beta_M[,1] = Beta_M[,1]
colnames(Shrinking_Beta_M) = colnames(Beta_M)

Shrinking_Beta_Up = data.frame(Shrinking_Beta_Up)
Shrinking_Beta_Up[,1] = Beta_Up[,1]
colnames(Shrinking_Beta_Up) = colnames(Beta_Up)

Shrinking_Beta_Down = data.frame(Shrinking_Beta_Down)
Shrinking_Beta_Down[,1] = Beta_Down[,1]
colnames(Shrinking_Beta_Down) = colnames(Beta_M)

Shrinking_Beta_CF = data.frame(Shrinking_Beta_CF)
Shrinking_Beta_CF[,1] = Beta_CF[,1]
colnames(Shrinking_Beta_CF) = colnames(Beta_CF)

Shrinking_Beta_DR = data.frame(Shrinking_Beta_DR)
Shrinking_Beta_DR[,1] = Beta_DR[,1]
colnames(Shrinking_Beta_DR) = colnames(Beta_DR)


load("Weight_Beta_M.RData")
load("Weight_Beta_Up.RData")
load("Weight_Beta_Down.RData")
load("Weight_Beta_CF.RData")
load("Weight_Beta_DR.RData")
#This is the Cross-section
Weighting_Shrinking_Beta_M = matrix(NA,nrow=dim(Weighting_Beta_M)[1],ncol=dim(Weighting_Beta_M)[2])
Weighting_Shrinking_Beta_Up = matrix(NA,nrow=dim(Weighting_Beta_M)[1],ncol=dim(Weighting_Beta_M)[2])
Weighting_Shrinking_Beta_Down = matrix(NA,nrow=dim(Weighting_Beta_M)[1],ncol=dim(Weighting_Beta_M)[2])
Weighting_Shrinking_Beta_CF = matrix(NA,nrow=dim(Weighting_Beta_M)[1],ncol=dim(Weighting_Beta_M)[2])
Weighting_Shrinking_Beta_DR = matrix(NA,nrow=dim(Weighting_Beta_M)[1],ncol=dim(Weighting_Beta_M)[2])
library(data.table)
setDF(Weighting_Beta_M)
setDF(Weighting_Beta_Up)
setDF(Weighting_Beta_Down)
setDF(Weighting_Beta_CF)
setDF(Weighting_Beta_DR)

#For Beta_M
cross_Beta_M = rowMeans(Weighting_Beta_M[,2:ncol(Weighting_Beta_M)], na.rm=T)
sd_cross_Beta_M = (sd(rowMeans(Weighting_Beta_M[,2:ncol(Weighting_Beta_M)], na.rm=T),na.rm=T))^2
prior_Beta_M = (apply(Weighting_Beta_M[,2:ncol(Weighting_Beta_M)],1,sd,na.rm=T))^2
#For Beta_Up
cross_Beta_Up = rowMeans(Weighting_Beta_Up[,2:ncol(Weighting_Beta_Up)], na.rm=T)
sd_cross_Beta_Up = (sd(rowMeans(Weighting_Beta_Up[,2:ncol(Weighting_Beta_Up)], na.rm=T),na.rm=T))^2
prior_Beta_Up = (apply(Weighting_Beta_Up[,2:ncol(Weighting_Beta_Up)],1,sd,na.rm=T))^2
#For Beta_Down
cross_Beta_Down = rowMeans(Weighting_Beta_Down[,2:ncol(Weighting_Beta_Down)], na.rm=T)
sd_cross_Beta_Down = (sd(rowMeans(Weighting_Beta_Down[,2:ncol(Weighting_Beta_Down)], na.rm=T),na.rm=T))^2
prior_Beta_Down = (apply(Weighting_Beta_Down[,2:ncol(Weighting_Beta_Down)],1,sd,na.rm=T))^2
#For Beta CF
cross_Beta_CF = rowMeans(Weighting_Beta_CF[,2:ncol(Weighting_Beta_CF)], na.rm=T)
sd_cross_Beta_CF = (sd(rowMeans(Weighting_Beta_CF[,2:ncol(Weighting_Beta_CF)], na.rm=T),na.rm=T))^2
prior_Beta_CF = (apply(Weighting_Beta_CF[,2:ncol(Weighting_Beta_CF)],1,sd,na.rm=T))^2
#For Beta DR
cross_Beta_DR = rowMeans(Weighting_Beta_DR[,2:ncol(Weighting_Beta_DR)], na.rm=T)
sd_cross_Beta_DR = (sd(rowMeans(Weighting_Beta_DR[,2:ncol(Weighting_Beta_DR)], na.rm=T),na.rm=T))^2
prior_Beta_DR = (apply(Weighting_Beta_DR[,2:ncol(Weighting_Beta_DR)],1,sd,na.rm=T))^2

#Do we shrink this based on the entire sample, apparently (...)
for(i in 2:ncol(Beta_M)){
  Weighting_Shrinking_Beta_M[,i] = (Weighting_Beta_M[,i] + sd_cross_Beta_M/prior_Beta_M * cross_Beta_M)/(1 + sd_cross_Beta_M/prior_Beta_M)
  
  Weighting_Shrinking_Beta_Up[,i] = (Weighting_Beta_Up[,i] + sd_cross_Beta_Up/prior_Beta_Up * cross_Beta_Up)/(1 + sd_cross_Beta_Up/prior_Beta_Up)
  
  Weighting_Shrinking_Beta_Down[,i] = (Weighting_Beta_Down[,i] + sd_cross_Beta_Down/prior_Beta_Down * cross_Beta_Down)/(1 + sd_cross_Beta_Down/prior_Beta_Down)
  
  Weighting_Shrinking_Beta_CF[,i] = (Weighting_Beta_CF[,i] + sd_cross_Beta_CF/prior_Beta_CF * cross_Beta_CF)/(1 + sd_cross_Beta_CF/prior_Beta_CF)
  
  Weighting_Shrinking_Beta_DR[,i] = (Weighting_Beta_DR[,i] + sd_cross_Beta_DR/prior_Beta_DR * cross_Beta_DR)/(1 + sd_cross_Beta_DR/prior_Beta_DR)
  
  print(paste0(i," of iteration ",ncol(Weighting_Beta_M)," which is ",i/ncol(Weighting_Beta_M)*100," Percent"))
}
#Add the Dates:
Weighting_Shrinking_Beta_M = data.frame(Weighting_Shrinking_Beta_M)
Weighting_Shrinking_Beta_M[,1] = Weighting_Beta_M[,1]
colnames(Weighting_Shrinking_Beta_M) = colnames(Weighting_Beta_M)

Weighting_Shrinking_Beta_Up = data.frame(Weighting_Shrinking_Beta_Up)
Weighting_Shrinking_Beta_Up[,1] = Weighting_Beta_Up[,1]
colnames(Weighting_Shrinking_Beta_Up) = colnames(Weighting_Beta_Up)

Weighting_Shrinking_Beta_Down = data.frame(Weighting_Shrinking_Beta_Down)
Weighting_Shrinking_Beta_Down[,1] = Weighting_Beta_Down[,1]
colnames(Weighting_Shrinking_Beta_Down) = colnames(Weighting_Beta_Down)

Weighting_Shrinking_Beta_CF = data.frame(Weighting_Shrinking_Beta_CF)
Weighting_Shrinking_Beta_CF[,1] = Weighting_Beta_CF[,1]
colnames(Weighting_Shrinking_Beta_CF) = colnames(Weighting_Beta_CF)

Weighting_Shrinking_Beta_DR = data.frame(Weighting_Shrinking_Beta_DR)
Weighting_Shrinking_Beta_DR[,1] = Weighting_Beta_DR[,1]
colnames(Weighting_Shrinking_Beta_DR) = colnames(Weighting_Beta_DR)



save(Shrinking_Beta_M,file="Shrinking_Beta_M.RData")
save(Shrinking_Beta_CF,file="Shrinking_Beta_CF.RData")
save(Shrinking_Beta_DR,file="Shrinking_Beta_DR.RData")
save(Shrinking_Beta_Up,file="Shrinking_Beta_Up.RData")
save(Shrinking_Beta_Down,file="Shrinking_Beta_Down.RData")


save(Weighting_Shrinking_Beta_M,file="Weighting_Shrinking_Beta_M.RData")
save(Weighting_Shrinking_Beta_CF,file="Weighting_Shrinking_Beta_CF.RData")
save(Weighting_Shrinking_Beta_DR,file="Weighting_Shrinking_Beta_DR.RData")
save(Weighting_Shrinking_Beta_Up,file="Weighting_Shrinking_Beta_Up.RData")
save(Weighting_Shrinking_Beta_Down,file="Weighting_Shrinking_Beta_Down.RData")

########################################################################
#############################################################################################
#################################### Naive Solution #######################################################
# with a considerably longer runtime. Running this code is necessary to get the IVOL.

Beta_CF = matrix(NA,nrow=dim(CrossSection[[3]])[1],ncol=dim(CrossSection[[3]])[2])
Beta_DR = matrix(NA,nrow=dim(CrossSection[[3]])[1],ncol=dim(CrossSection[[3]])[2])
Beta_M  = matrix(NA,nrow=dim(CrossSection[[3]])[1],ncol=dim(CrossSection[[3]])[2])
#Check with Beta_M2
Beta_M2 = matrix(NA,nrow=dim(CrossSection[[3]])[1],ncol=dim(CrossSection[[3]])[2])
#Up and Downside Beta
Beta_Down = matrix(NA,nrow=dim(CrossSection[[3]])[1],ncol=dim(CrossSection[[3]])[2])
Beta_Up   = matrix(NA,nrow=dim(CrossSection[[3]])[1],ncol=dim(CrossSection[[3]])[2])

Weight_Beta_CF = matrix(NA,nrow=dim(CrossSection[[3]])[1],ncol=dim(CrossSection[[3]])[2])
Weight_Beta_DR = matrix(NA,nrow=dim(CrossSection[[3]])[1],ncol=dim(CrossSection[[3]])[2])
Weight_Beta_M  = matrix(NA,nrow=dim(CrossSection[[3]])[1],ncol=dim(CrossSection[[3]])[2])
#Check with Beta_M2
Weight_Beta_M2 = matrix(NA,nrow=dim(CrossSection[[3]])[1],ncol=dim(CrossSection[[3]])[2])
#Up and Downside Beta
Weight_Beta_Down = matrix(NA,nrow=dim(CrossSection[[3]])[1],ncol=dim(CrossSection[[3]])[2])
Weight_Beta_Up   = matrix(NA,nrow=dim(CrossSection[[3]])[1],ncol=dim(CrossSection[[3]])[2])


#IVOL
I_Vol = matrix(NA,nrow=dim(CrossSection[[3]])[1],ncol=dim(CrossSection[[3]])[2])




#Rolling Window of 5 years standard FM (1973), see Hollstein et al. (2019)
options(warn=2)
for(i in 2:dim(CrossSection[[3]])[2]){
  for(j in 61:dim(CrossSection[[3]])[1]){
    #Weighting:
      Weight_Beta_CF[j,i] = cov(Weights*(CrossSection[[3]][(j-59):j,i]-RF[(j-59):j]),Weights*N_CF[(j-59):j],use="pairwise.complete.obs")/var((Weights*(N_CF[(j-59):j]-N_DR[(j-59):j]))[!is.na(CrossSection[[3]][(j-59):j,i])],na.rm=T)+
        cov(Weights*(CrossSection[[3]][(j-59):j,i]-RF[(j-59):j]),Weights*N_CF[(j-60):(j-1)],use="pairwise.complete.obs")/var((Weights*(N_CF[(j-59):j]-N_DR[(j-59):j]))[!is.na(CrossSection[[3]][(j-59):j,i])],na.rm=T)
      
      Weight_Beta_DR[j,i] = cov(Weights*(CrossSection[[3]][(j-59):j,i]-RF[(j-59):j]),Weights*N_DR[(j-59):j],use="pairwise.complete.obs")/var((Weights*(N_CF[(j-59):j]-N_DR[(j-59):j]))[!is.na(CrossSection[[3]][(j-59):j,i])],na.rm=T)+
        cov(Weights*(CrossSection[[3]][(j-59):j,i]-RF[(j-59):j]),Weights*N_DR[(j-60):(j-1)],use="pairwise.complete.obs")/var((Weights*(N_CF[(j-59):j]-N_DR[(j-59):j]))[!is.na(CrossSection[[3]][(j-59):j,i])],na.rm=T)
      
      Weight_Beta_M[j,i] = Beta_CF[j,i] + Beta_DR[j,i]
      Weight_Beta_M2[j,i] = cov(Weights*(CrossSection[[3]][(j-59):j,i]-RF[(j-59):j]),Weights*CheckCheck$R_Me[(j-59):j],use="pairwise.complete.obs")/var(Weights[!is.na(CrossSection[[3]][(j-59):j,i]-RF[(j-59):j])]*(CheckCheck$R_Me[(j-59):j]-RF[(j-59):j])[!is.na(CrossSection[[3]][(j-59):j,i]-RF[(j-59):j])],na.rm=T)  
    
      #Beta Downside Ang et al. (2006) - Downside Risk
      if(length(CrossSection[[3]][(j-59):j,i][CheckCheck$R_Me[(j-59):j]<RF[(j-59):j]])!=0){
        Weight_Beta_Down[j,i] = cov(Weights[CheckCheck$R_Me[(j-59):j]<RF[(j-59):j]]*(CrossSection[[3]][(j-59):j,i][CheckCheck$R_Me[(j-59):j]<RF[(j-59):j]]-RF[(j-59):j][CheckCheck$R_Me[(j-59):j]<RF[(j-59):j]]),Weights[CheckCheck$R_Me[(j-59):j]<RF[(j-59):j]]*CheckCheck$R_Me[(j-59):j][CheckCheck$R_Me[(j-59):j]<RF[(j-59):j]],use="pairwise.complete.obs")/var((Weights[CheckCheck$R_Me[(j-59):j]<RF[(j-59):j]]*CheckCheck$R_Me[(j-59):j][CheckCheck$R_Me[(j-59):j]<RF[(j-59):j]])[!is.na(CrossSection[[3]][(j-59):j,i][CheckCheck$R_Me[(j-59):j]<RF[(j-59):j]])],na.rm=T)
      }else{
        Weight_Beta_Down[j,i] = NA
      }
      if(length(CrossSection[[3]][(j-59):j,i][CheckCheck$R_Me[(j-59):j]>RF[(j-59):j]])!=0){
      #Beta Upside Ang et al. (2006) - Downside Risk
        Weight_Beta_Up[j,i] = cov(Weights[CheckCheck$R_Me[(j-59):j]>RF[(j-59):j]]*(CrossSection[[3]][(j-59):j,i][CheckCheck$R_Me[(j-59):j]>RF[(j-59):j]]-RF[(j-59):j][CheckCheck$R_Me[(j-59):j]>RF[(j-59):j]]),Weights[CheckCheck$R_Me[(j-59):j]>RF[(j-59):j]]*CheckCheck$R_Me[(j-59):j][CheckCheck$R_Me[(j-59):j]>RF[(j-59):j]],use="pairwise.complete.obs")/var((Weights[CheckCheck$R_Me[(j-59):j]>RF[(j-59):j]]*CheckCheck$R_Me[(j-59):j][CheckCheck$R_Me[(j-59):j]>RF[(j-59):j]])[!is.na(CrossSection[[3]][(j-59):j,i][CheckCheck$R_Me[(j-59):j]>RF[(j-59):j]])],na.rm=T)
      }else{
        Weight_Beta_Up[j,i] = NA
      }
      #Historical Beta
      #With the news series, we can calculate the betas easily, by:
      Beta_CF[j,i] = cov((CrossSection[[3]][(j-59):j,i]-RF[(j-59):j]),N_CF[(j-59):j],use="pairwise.complete.obs")/var((N_CF[(j-59):j] - N_DR[(j-59):j])[!is.na(CrossSection[[3]][(j-59):j,i]-RF[(j-59):j])],na.rm=T)+
        cov((CrossSection[[3]][(j-59):j,i]-RF[(j-59):j]),N_CF[(j-60):(j-1)],use="pairwise.complete.obs")/var((N_CF[(j-59):j] - N_DR[(j-59):j])[!is.na(CrossSection[[3]][(j-59):j,i]-RF[(j-59):j])],na.rm=T)
    
      Beta_DR[j,i] = cov((CrossSection[[3]][(j-59):j,i]-RF[(j-59):j]),-N_DR[(j-59):j],use="pairwise.complete.obs")/var((N_CF[(j-59):j] - N_DR[(j-59):j])[!is.na(CrossSection[[3]][(j-59):j,i]-RF[(j-59):j])],na.rm=T)+
        cov((CrossSection[[3]][(j-59):j,i]-RF[(j-59):j]),-N_DR[(j-60):(j-1)],use="pairwise.complete.obs")/var((N_CF[(j-59):j] - N_DR[(j-59):j])[!is.na(CrossSection[[3]][(j-59):j,i]-RF[(j-59):j])],na.rm=T)
    
      #"Regular" Beta is then simply the sum of these two, see Campbell and Vuolteenaho (2004):
      Beta_M[j,i] = Beta_CF[j,i] + Beta_DR[j,i]
      #Check with:
      Beta_M2[j,i] = cov(CrossSection[[3]][(j-59):j,i]-RF[(j-59):j],CheckCheck$R_Me[(j-59):j],use="pairwise.complete.obs")/var((CheckCheck$R_Me[(j-59):j])[!is.na(CrossSection[[3]][(j-59):j,i])],na.rm=T)
    
      if(length(CrossSection[[3]][(j-59):j,i][CheckCheck$R_Me[(j-59):j]<RF[(j-59):j]])!=0){
      #Beta Downside Ang et al. (2006) - Downside Risk
        Beta_Down[j,i] = cov(CrossSection[[3]][(j-59):j,i][CheckCheck$R_Me[(j-59):j]<RF[(j-59):j]]-RF[(j-59):j][CheckCheck$R_Me[(j-59):j]<RF[(j-59):j]],CheckCheck$R_Me[(j-59):j][CheckCheck$R_Me[(j-59):j]<RF[(j-59):j]],use="pairwise.complete.obs")/var((CheckCheck$R_Me[(j-59):j][CheckCheck$R_Me[(j-59):j]<RF[(j-59):j]])[!is.na(CrossSection[[3]][(j-59):j,i][CheckCheck$R_Me[(j-59):j]<RF[(j-59):j]])],na.rm=T)
      }else{
        Beta_Down[j,i] = NA
      }
      if(length(CrossSection[[3]][(j-59):j,i][CheckCheck$R_Me[(j-59):j]>RF[(j-59):j]])!=0){
      #Beta Upside Ang et al. (2006) - Downside Risk
        Beta_Up[j,i] = cov(CrossSection[[3]][(j-59):j,i][CheckCheck$R_Me[(j-59):j]>RF[(j-59):j]]-RF[(j-59):j][CheckCheck$R_Me[(j-59):j]>RF[(j-59):j]],CheckCheck$R_Me[(j-59):j][CheckCheck$R_Me[(j-59):j]>RF[(j-59):j]],use="pairwise.complete.obs")/var((CheckCheck$R_Me[(j-59):j][CheckCheck$R_Me[(j-59):j]>RF[(j-59):j]])[!is.na(CrossSection[[3]][(j-59):j,i][CheckCheck$R_Me[(j-59):j]>RF[(j-59):j]])],na.rm=T)
      }else{
        Beta_Up[j,i] = NA
      }
    if(sum(is.na((CrossSection[[3]][(j-59):j,i]-RF[(j-59):j])))<6){
      #Calculate the residuals for the return series:
      library(pracma)
      Temp = na.omit(cbind(CrossSection[[3]][(j-59):j,i]-RF[(j-59):j],rep(1,60),FF[(j-59):j,c("Mkt.RF","SMB","HML")]))
      coeff_ret = mldivide(as.matrix(Temp[,2:5]),as.matrix(Temp[,1]))
      #Calculate the residuals:
      data_res1 = (Temp[,1]) - as.matrix(Temp[,2:5])%*%coeff_ret
      I_Vol[j,i] = sd(data_res1,na.rm=T)
    }
    }
  if(i%%1000 == 0){
    print(paste("iteration", i, "of", dim(CrossSection[[3]])[2]))
  }
}
options(warn=1)

#Estimate With Shrinkage:
#i Companies
#j times

Shrinking_Beta_CF = matrix(NA,nrow=dim(CrossSection[[3]])[1],ncol=dim(CrossSection[[3]])[2])
Shrinking_Beta_DR = matrix(NA,nrow=dim(CrossSection[[3]])[1],ncol=dim(CrossSection[[3]])[2])
Shrinking_Beta_M  = matrix(NA,nrow=dim(CrossSection[[3]])[1],ncol=dim(CrossSection[[3]])[2])
#Check with Beta_M2
Shrinking_Beta_M2 = matrix(NA,nrow=dim(CrossSection[[3]])[1],ncol=dim(CrossSection[[3]])[2])
#Up and Downside Beta
Shrinking_Beta_Down = matrix(NA,nrow=dim(CrossSection[[3]])[1],ncol=dim(CrossSection[[3]])[2])
Shrinking_Beta_Up   = matrix(NA,nrow=dim(CrossSection[[3]])[1],ncol=dim(CrossSection[[3]])[2])


std <- function(x, ...) sd(x, ...)/sqrt(length(x))

for(i in 2:dim(CrossSection[[3]])[2]){
  for(j in 12:dim(CrossSection[[3]])[1]){
    if(sum(is.na(Beta_CF[(j-11):j,i]))<6){
      #Shrinking:
      Shrinking_Beta_CF[j,i]   = (Beta_CF[j,i] + ( std(Beta_CF[(j-11):j,i],na.rm=T)^2/std(Beta_CF[j,2:dim(CrossSection[[3]])[2]],na.rm=T)^2 * mean(Beta_CF[j,2:dim(CrossSection[[3]])[2]],na.rm=T)))/
  (1 + ( std(Beta_CF[(j-11):j,i],na.rm=T)^2/std(Beta_CF[j,2:dim(CrossSection[[3]])[2]],na.rm=T)^2))
    }else{
      Shrinking_Beta_CF[j,i]   = NA
}
    if(sum(is.na(Beta_DR[(j-11):j,i]))<6){
Shrinking_Beta_DR[j,i]   = (Beta_DR[j,i] + ( std(Beta_DR[(j-11):j,i],na.rm=T)^2/std(Beta_DR[j,2:dim(CrossSection[[3]])[2]],na.rm=T)^2 * mean(Beta_DR[j,2:dim(CrossSection[[3]])[2]],na.rm=T)))/
  (1 + ( std(Beta_DR[(j-11):j,i],na.rm=T)^2/std(Beta_DR[j,2:dim(CrossSection[[3]])[2]],na.rm=T)^2))
    }else{
      Shrinking_Beta_DR[j,i]   = NA
}
    if(sum(is.na(Beta_M[(j-11):j,i]))<6){
Shrinking_Beta_M[j,i]    = (Beta_M[j,i] + ( std(Beta_M[(j-11):j,i],na.rm=T)^2/std(Beta_M[j,2:dim(CrossSection[[3]])[2]],na.rm=T)^2 * mean(Beta_M[j,2:dim(CrossSection[[3]])[2]],na.rm=T)))/
  (1 + ( std(Beta_M[(j-11):j,i],na.rm=T)^2/std(Beta_M[j,2:dim(CrossSection[[3]])[2]],na.rm=T)^2))
    }else{
      Shrinking_Beta_M[j,i]    = NA
    }
    if(sum(is.na(Beta_M2[(j-11):j,i]))<6){
Shrinking_Beta_M2[j,i]   = (Beta_M2[j,i] + ( std(Beta_M2[(j-11):j,i],na.rm=T)^2/std(Beta_M2[j,2:dim(CrossSection[[3]])[2]],na.rm=T)^2 * mean(Beta_M2[j,2:dim(CrossSection[[3]])[2]],na.rm=T)))/
  (1 + ( std(Beta_M2[(j-11):j,i],na.rm=T)^2/std(Beta_M2[j,2:dim(CrossSection[[3]])[2]],na.rm=T)^2))
    }else{
      Shrinking_Beta_M2[j,i]   = NA 
    }
    if(sum(is.na(Beta_Down[(j-11):j,i]))<6){
Shrinking_Beta_Down[j,i] = (Beta_Down[j,i] + ( std(Beta_Down[(j-11):j,i],na.rm=T)^2/std(Beta_Down[j,2:dim(CrossSection[[3]])[2]],na.rm=T)^2 * mean(Beta_Down[j,2:dim(CrossSection[[3]])[2]],na.rm=T)))/
  (1 + ( std(Beta_Down[(j-11):j,i],na.rm=T)^2/std(Beta_Down[j,2:dim(CrossSection[[3]])[2]],na.rm=T)^2))
    }else{
      Shrinking_Beta_Down[j,i] = NA
    }
    if(sum(is.na(Beta_Up[(j-11):j,i]))<6){
Shrinking_Beta_Up[j,i]   = (Beta_Up[j,i] + ( std(Beta_Up[(j-11):j,i],na.rm=T)^2/std(Beta_Up[j,2:dim(CrossSection[[3]])[2]],na.rm=T)^2 * mean(Beta_Up[j,2:dim(CrossSection[[3]])[2]],na.rm=T)))/
  (1 + ( std(Beta_Up[(j-11):j,i],na.rm=T)^2/std(Beta_Up[j,2:dim(CrossSection[[3]])[2]],na.rm=T)^2))
    }else{
      Shrinking_Beta_Up[j,i]   = NA
}
  }
  if(i%%1000 == 0){
    print(paste("iteration", i, "of", dim(CrossSection[[3]])[2]))
  }
}


which(colSums(!is.na(Beta_M)) > 400)

mean(Beta_M,na.rm=T)
quantile(Beta_M, probs  = c(0.000001,0.001,0.01,0.05,0.25,0.5,0.75,0.95,0.99), na.rm=T)
mean(Beta_M2,na.rm=T)
quantile(Beta_M2, probs = c(0.000001,0.001,0.01,0.05,0.25,0.5,0.75,0.95,0.99), na.rm=T)
mean(Beta_CF,na.rm=T)
quantile(Beta_CF, probs = c(0.000001,0.001,0.01,0.05,0.25,0.5,0.75,0.95,0.99), na.rm=T)
mean(Beta_DR,na.rm=T)
quantile(Beta_DR, probs = c(0.000001,0.001,0.01,0.05,0.25,0.5,0.75,0.95,0.99), na.rm=T)
mean(Beta_Up,na.rm=T)
quantile(Beta_Up, probs = c(0.000001,0.001,0.01,0.05,0.25,0.5,0.75,0.95,0.99), na.rm=T)
mean(Beta_Down,na.rm=T)
quantile(Beta_Down, probs = c(0.000001,0.001,0.01,0.05,0.25,0.5,0.75,0.95,0.99), na.rm=T)


################################################ Save the results ##############################################################




#Add the Dates from CrossSection[[3]] to Beta_CF, Beta_DR and Beta_M
Beta_M[,1]  = CrossSection[[3]][,1]
Beta_M2[,1] = CrossSection[[3]][,1]
Beta_CF[,1] = CrossSection[[3]][,1]
Beta_DR[,1] = CrossSection[[3]][,1]
Beta_Up[,1] = CrossSection[[3]][,1]
Beta_Down[,1] = CrossSection[[3]][,1]
Weight_Beta_M[,1]  = CrossSection[[3]][,1]
Weight_Beta_M2[,1] = CrossSection[[3]][,1]
Weight_Beta_CF[,1] = CrossSection[[3]][,1]
Weight_Beta_DR[,1] = CrossSection[[3]][,1]
Weight_Beta_Up[,1] = CrossSection[[3]][,1]
Weight_Beta_Down[,1] = CrossSection[[3]][,1]
Shrinking_Beta_M[,1] = CrossSection[[3]][,1]
Shrinking_Beta_M2[,1] = CrossSection[[3]][,1]
Shrinking_Beta_CF[,1] = CrossSection[[3]][,1]
Shrinking_Beta_DR[,1] = CrossSection[[3]][,1]
Shrinking_Beta_Up[,1] = CrossSection[[3]][,1]
Shrinking_Beta_Down[,1] = CrossSection[[3]][,1]
I_Vol[,1] = CrossSection[[3]][,1]
#Add Sec ID's
Beta_M  = rbind(c(Sec_ID),Beta_M)
Beta_M2 = rbind(c(Sec_ID),Beta_M2)
Beta_CF = rbind(c(Sec_ID),Beta_CF)
Beta_DR = rbind(c(Sec_ID),Beta_DR)
Beta_Up = rbind(c(Sec_ID),Beta_Up)
Beta_Down = rbind(c(Sec_ID),Beta_Down)
Weight_Beta_M  = rbind(c(Sec_ID),Weight_Beta_M)
Weight_Beta_M2 = rbind(c(Sec_ID),Weight_Beta_M2)
Weight_Beta_CF = rbind(c(Sec_ID),Weight_Beta_CF)
Weight_Beta_DR = rbind(c(Sec_ID),Weight_Beta_DR)
Weight_Beta_Up = rbind(c(Sec_ID),Weight_Beta_Up)
Weight_Beta_Down = rbind(c(Sec_ID),Weight_Beta_Down)
Shrinking_Beta_M  = rbind(c(Sec_ID),Shrinking_Beta_M)
Shrinking_Beta_M2 = rbind(c(Sec_ID),Shrinking_Beta_M2)
Shrinking_Beta_CF = rbind(c(Sec_ID),Shrinking_Beta_CF)
Shrinking_Beta_DR = rbind(c(Sec_ID),Shrinking_Beta_DR)
Shrinking_Beta_Up = rbind(c(Sec_ID),Shrinking_Beta_Up)
Shrinking_Beta_Down = rbind(c(Sec_ID),Shrinking_Beta_Down)
I_Vol = rbind(c(Sec_ID),I_Vol)

save(Beta_M,file="Beta_M.RData")
save(Beta_M2,file="Beta_M2.RData")
save(Beta_CF,file="Beta_CF.RData")
save(Beta_DR,file="Beta_DR.RData")
save(Beta_Up,file="Beta_Up.RData")
save(Beta_Down,file="Beta_Down.RData")
save(Weight_Beta_M,file="Weight_Beta_M.RData")
save(Weight_Beta_M2,file="Weight_Beta_M2.RData")
save(Weight_Beta_CF,file="Weight_Beta_CF.RData")
save(Weight_Beta_DR,file="Weight_Beta_DR.RData")
save(Weight_Beta_Up,file="Weight_Beta_Up.RData")
save(Weight_Beta_Down,file="Weight_Beta_Down.RData")
save(Shrinking_Beta_M,file="Shrinking_Beta_M.RData")
save(Shrinking_Beta_M2,file="Shrinking_Beta_M2.RData")
save(Shrinking_Beta_CF,file="Shrinking_Beta_CF.RData")
save(Shrinking_Beta_DR,file="Shrinking_Beta_DR.RData")
save(Shrinking_Beta_Up,file="Shrinking_Beta_Up.RData")
save(Shrinking_Beta_Down,file="Shrinking_Beta_Down.RData")
save(I_Vol,file="I_Vol.RData")